install.packages("Seurat")
library(Seurat)
library(dplyr)
pbmc.data<-Read10X(data.dir="E:/R Studio/shixi/shixi 7/pbmc3k_filtered_gene_bc_matrices.tar/pbmc3k_/filtered_gene_bc_matrices/hg19/") #写入文件
pbmc<-CreateSeuratObject(counts=pbmc.data,project="pbmc3k",min.cells = 3,min.features = 200) #保留至少3个细胞中检测到表达基因，单个细胞中至少检测到200个基因
pbmc[['percent.mt']]<-PercentageFeatureSet(pbmc,pattern = "^MT-|^mt-") #人类用MT，小鼠用mt
VlnPlot(pbmc,features=c("nFeature_RNA","percent.mt"),ncol=3) #画出每个单细胞中的基因数量（nFeature_RNA）和线粒体基因百分比（percent.mt）的小提琴图
pbmc<-subset(pbmc,subset = nFeature_RNA>200&nFeature_RNA<2500&percent.mt<5)#从数据集中移除不需要的细胞，单个细胞中至少检测到200个基因且不多于2500个基因，线粒体基因比例不超过5%
pbmc<-NormalizeData(pbmc,normalization.method="LogNormalize",scale.fator=10000)#标准化数据，使用“LogNormalize”全局缩放的归一化方法
all.genes<-rownames(pbmc) #将pbmc矩阵的行名（即基因的名称）提取出来存储在变量all.genes中
pbmc<-ScaleData(pbmc,features=all.genes) #将所有基因执行数据缩放（scaling）操作并将结果保存回 pbmc 中
pbmc<-FindVariableFeatures(pbmc,selection.method="vst",nfeatures=2000)#识别2000个高变异基因并将这些基因的信息保存在 pbmc 对象中
pbmc<-ScaleData(pbmc) #对pbmc进行缩放操作
top10<-head(VariableFeatures(pbmc),10)#前10的高变异基因
plot1<-VariableFeaturePlot(pbmc)#绘制高变异基因分布图
plot2<-LabelPoints(plot=plot1,points=top10,repel=TRUE, xnudge=0, ynudge=0)#在高变异基因分布图中将前top10的高变异基因标注出来。xnudge=0, ynudge=0这个是系统建议加上去的
plot1 #可视化图表
plot2
pbmc<-RunPCA(pbmc,features = VariableFeatures(object = pbmc))#运行PCA主成分分析线性降维,指定了使用高变异基因
PCAPlot(pbmc) #主成分分析图表
PCA_plot<-DimPlot(pbmc,dims=1:2)#使用DimPlot函数可视化PCA结果。选择前2个维度进行查看dims=2指定了前两个主成分PC1和PC2进行可视化 此源于 https://www.jianshu.com/p/75c6d55a63bb
PCA_plot
print(pbmc[["pca"]],dims=1:5,nfeatures = 5)#选择前5个维度进行查看 此源于 https://www.jianshu.com/p/75c6d55a63bb
DimPlot(pbmc,reduction="pca")#可视化
pbmc<-FindNeighbors(pbmc,dims=1:10)#计算最邻近距离
pbmc<-FindClusters(pbmc,resolution=0.5)#进行聚类，设置了下游聚类的“间隔尺度”的分辨率参数resolution未0.5
head(Idents(pbmc),5)#查看前5个聚类细胞的id
pbmc<-RunUMAP(pbmc,dims=1:10)#进行UMAP非线性降维
UMAPPlot(pbmc)#进行可视化
pbmc.markers<-FindAllMarkers(pbmc,only.pos=TRUE,min.pct=0.25,logfc.threshold = 0.25)#找到每个cluster相比于其余cluster的markgene，只报道阳性的markgene
top_markers<-pbmc.markers%>%  #提取每个cluster的Top 2基因的名字
  group_by(cluster)%>%
  slice_max(n=2,order_by=avg_log2FC)%>%
  ungroup()%>%
  pull(gene)%>%
  unique()  #如果存在重复的基因名，使用unique()去除  
DotPlot(pbmc, features = top_markers) + RotatedAxis() #画出气泡图
VlnPlot(pbmc,features = c("NPM1","RPS27","CD79A"))#画出基因MS4A1、RPS27和CD79A的小提琴图
FeaturePlot(pbmc,features = c("MS4A1","RPS27","RPL9","LDHB")) #画出marker基因MS4A1、RPS27、RPL9和LDHB的feature plot
